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Axisymmetric  Analysis  of  the  Cylindrically 
Orthotropic  Disk  of  Variable  Fiber  Orientation 

Abstract 

Injection  molding  of  axisymmetric  bodies  with 
fiber-reinforced  molding  compounds  results  in  cylindrically 
orthotropic  components  in  which  the  fiber  orientation  varies 
with  radial  position.  Consequently,  development  of  analysis 
techniques  for  determining  the  effect  of  material  property 
variability  on  the  response  of  these  components  is  of 
great  importance.  In  the  present  study,  a  numerical 
integration  scheme  for  the  analysis  of  cylindrically 
orthotropic  annular  disks  with  variable  elastic  constants 
is  presented.  The  influence  of  material  property  variations 
temperature  changes  and  centrifugal  forces  on  the  response 
of  an  annular  disk  subjected  to  internal  and  external 
pressure  or  displacement  boundary  condition  is  included  in 
the  analysis i  Correlation  of  numerical  integration  results 
with  analytical  and  finite-element  solutions  is  excellent. 


Axisymmetric  Analysis  of  the  Cylindrically 

Orthotropic  Disk  of  Variable  Fiber  Orientation 

Introduction 

Injection  molding  of  axisymmetric  bodies  with  fiber- 
reinforced  molding  compounds  results  in  cylindrically 
orthotropic  components.  In  general,  the  fiber  orientation 
will  not  be  constant  throughout  the  body,  but  will  be  a 
function  of  radial  position.  It  has  recently  been  shown 
(i]  that  the  flow  field  and  mold  geometry  determine  fiber 
orientation.  For  example,  the  influence  of  converging 
and  diverging  flow  fields  on  fiber  orientation  are  shown 
clearly  in  Figure  1.  Consequently,  injection  molding  of 
an  axisymmetric  mold  uniformly  on  the  inner  or  outer  radius 
(axisymmetric  flow)  will  result  in  fiber  orientation  which 
is  dependent  solely  on  the  radial  position.  Material 
properties  are  determined  uniquely  by  fiber  orientation, 
constituent  properties,  fiber  aspect  ratio  and  fiber  volume 
fraction.  McCullough  [1,  2]  has  introduced  the  following 
parameters  to  quantify  fiber  orientation: 

f  =  1/2  [3  <  cos^(|)  >  -  1] 

g  =  1/4  [5  <  cos^(|)  >  -  1] 

Tr/2 

<cos  (|)  >  =  /o  N  ((())  cos^(|)sin(j)d({) 

where  (|)  is  the  angle  which  a  fiber  makes  with  the  longitudinal 


direction  and  N((j))  is  the  percentage  of  fibers  with  that 
direction.  In  figure  2,  the  bounds  on  Young's  modulus 
as  a  function  of  the  orientation  parameter  "f"  are  shown. 
Clearly,  the  variation  of  fiber  orientation  with  radial 
position  significantly  influences  the  mechanical  properties 
of  the  component. 

In  the  present  study,  a  numerical  integration  scheme 
is  developed  to  determine  the  response  of  cylindrically 
orthotropic  disks  with  elastic  constants  which  vary  in  the 
radial  direction.  Uniform  temperature  variations  and 
centrifugal  body  forces,  as  well  as,  pressure  on  displacement 
boundary  conditions  prescribed  on  the  inner  and  outer  radii 
are  considered. 

Formulation  of  the  governing  equation  in  terms  of 
displacements  yields  a  second  order  linear  ordinary  differen¬ 
tial  equation  with  non-constant  coefficients  (See  Appendix  A) . 
In  general,  a  closed  form  solution  of  this  equation  does  not 
exist.  The  integration  scheme  requires  the  reduction  of 
the  governing  equation  to  two  simultaneous  first  order 
differential  equations  which  are  solved  using  Hamming's 
Predictor-Corrector  Method  [3].  Unfortunately,  Hamming's 
Method  requires  two  initial  value  conditions  (one  for  each 
first  order  equation)  whereas  only  one  is  known  in  the  actual 
boundary  value  problem.  Consequently,  a  half-interval  search 
technique  is  incorporated  into  the  program  in  which  upper 
and  lower  bounds  for  the  unknown  initial  condition  are 


prescribed.  The  average  value  is  employed  in  the  integration 
and  correlation  of  the  solution  with  the  second  known  boundary 
condition  enables  the  interval  of  uncertainty  to  be  halved. 
Subsequent  iterations  converge  quickly  to  the  solution. 

In  fact,  if  is  the  length  of  the  starting  interval,  then 
the  number  (N)  of  interval  halving  operations  required  to 
reduce  the  interval  of  uncertainty  to  Aj^  is  given  by 


Closed-form  solutions  for  two  special  variations 
of  elastic  properties  are  presented  in  Appendices  B  and  C 
to  verify  numerical  integration  results.  The  first  solution 
is  for  uniform  properties  and  the  second  assumes  that  the 
modulii  vary  along  a  radius  according  to  a  power  law; 


E  = 
r 


= 


where  m  is  an  arbitrary  real  number,  and  the  Poisson's 
ratios  are  held  constant.  In  addition,  finite  element 
results  for  a  linear  variation  of  properties  are  compared 
to  the  numerical  integration  results . 


Correlation  of  Results 


Agreement  of  numerical  integration  results  with 
analytical  and  finite  element  results  was  found  to  be 
excellent.  In  Figures  3-6,  typical  results  for  temperature 
variations,  centrifugal  body  forces  and  internal  and  external 
radial  stress  tractions  are  presented.  The  solid  line  in 
all  figures  corresponds  to  the  numerical  integration  prediction. 
Superimposed  are  analytic  and  finite  element  results  shown 
as  symbols.  The  excellent  agreement  is  obtained  using 
an  integration  step-size  of  0.001  inches  (0.025  mm). 
Approximately  15-20  iterations  are  required  to  determine  the 
unknown  initial  condition  with  sufficient  accuracy  to 
satisfy  the  remaining  boundary  condition  on  the  outer  radius. 
Although  only  the  correlation  of  stress  components  are 
presented  in  Figures  3-6,  excellent  agreement  was 
obtained  for  displacement  and  strain  component  values  as  well. 

For  reference,  the  material  properties  employed  in 
the  various  solutions  are  given  below: 

Constant  Properties:  Radial  Fiber  Orientation 

=  2.6  Msi(17.9GPa)  =  4.5xlO"®/°F  ( 8 . Ixl0"®/°C) 

Eg  =  1.4  Msi(9.65GPa)  aQ=  10.1xl0”®/°F  ( 18 . 2xl0" V°C) 

%r  =  0-27 

Constant  Properties:  Tangential  Fiber  Orientation 

E^  =  1.4  Msi(9.65GPa)  =  10.1xl0“^/°F  ( 18 . 2xlO~®/°C) 

E  g  =  2.6Mai  (17.9GPa)  ag=  4.5xlO"®/°F  ( 8  .  Ixl0“^/°C) 

er 


V 


.501 


Power  Law  Variation  (m  specified) 


=  2.6  r™  Msi  (IT.Sr'^GPa)  v^q  =  .501 
E.  =  1.4  r”*  Msi  (9.65r™GPa)  a  =  (4 . 5xl0“®)  r”‘/°F 

U  IT 

(8.1xlO“®r"'/°C) 

Vqj.  =  .27  Wq  =  (10.1x10“^)  r"‘/°F 

(18.2xlO"®r"‘/°C) 


Linear  Variation 

Er  =  (r-a)  +  2.6)Msi  ( (— (r-a)  +  17.9)  GPa) 

Eq  =  (r-a)  +  1.4) Msi  ( (|^)  (r-a)  +  9.65)  GPa) 


^er  = 


^re  ~  '^0r  ®r/^0 


=  +  4.5)xlO"^/°F  ( (— tl-IIT"-  +  8.1)xl0"^/ 


(b-a) 


“0  ""  (-■^(b-a?^^  lO.DxlO  ^/°F  (  +  I8.2)xl0  V°C) 


where 

E^  -  Young's  Modulus  in  radial  direction 
Eg  -  Young's  Modulus  in  tangential  direction 
Vr0  Vg^  -  Poisson  ratios 

-  Thermal  coefficient  of  expansion  in  radial  direction 
Ug'-  Thermal  coefficient  of  expansion  in  tangential  direction 


User's  Guide 


Fortran  computer  codes  have  been  developed  for  two 
analytical  solutions  and  for  the  Hamming's  Predictor- 
Corrector  Method.  Analysis  details  and  program  listings 
may  be  found  in  Appendices  A,  B  and  C.  The  programs  have 
been  written  in  an  interactive  format  which  necessitates 
execution  from  a  terminal  or  similar  device.  In  Table  1, 
program  symbol  definitions  are  defined.  The  following 
examples  will  be  illustrative. 


Analytic  Solution;  Constant  Properties  (VARPROP/CFD) 

Table  2  indicates  the  line  numbers  in  VARPROP/CFD 
which  describe  the  material  properties .  These  are  the 
only  lines  that  must  be  altered  when  another  set  of  properties 
are  to  be  input.  In  Figure  7  a  sample  program  executioh  is 
presented  where  data  input  is  requested  by  the  program. 

Note  that  displacement  boundary  conditions  are  not  possible. 


Analytic  Solution;  Power  Law  Variation  (VARPR0P/CF2) 

The  anlytic  solution  assumes  the  following  material 


property  variation: 

i-i  ro 

=  E  r 

r  rm 

a  =  a  r”^ 
r  rm 

Vg^  =  constant 

m 

“B  =  “Bm^ 

In  Table  3,  the  line  numbers  in  VARPR0P/CF2  which  describe 
the  material  properties  are  shown.  Execution  is  straight¬ 
forward  as  indicated  in  Figure  8.  Note  that  centrifugal 
body  forces  and  displacement  boundary  conditions  are  not 
included. 

Numerical  Integration  (VARPROP/NUMD) 

The  material  property  section  of  VARPROP/WUMD  is 
located  in  Subroutine  PROP  (See  Appendix  A  for  further 
detail) .  All  property  dependence  on  radial  position  must 
be  input.  Derivatives  of  these  properties  are  also  required. 

In  Table  4,  for  example,  the  material  properties  employed 
in  the  correlation  of  numerical  with  analytical  results 
(constant  properties)  are  presented.  Material  property 
input  for  power  law  and  linear  variations  are  shown  in 
Tables  5  and  6,  respectively.  Note  that  derivatives  are 
input  in  a  straightforward  manner. 

Execution  of  VARPROP/NuMD  is  also  in  an  interactive 
mode.  The  program  versatility  allows  displacement  o^ 
stress  boundary  conditions  in  the  inner  and  outer  surface. 

As  mentioned  previously,  upper  and  lower  bounds  on  the 
unknown  initial  condition  at  the  inner  radius  must  be 
supplied.  If  a  radial  stress  initial  condition  is  specified, 
bounds  on  the  tangential  stress  components  are  input  as 
shown  in  Figures  9  and  10.  For  displacement  initial  conditions, 
bounds  on  the  radial  strain  component  at  the  inner  radius 


are  required  (See  Figures  11  and  12) .  Results  from  numerical 
integration  and  the  analytical  solution  for  constant  properties 
(radial  fiber  orientation)  are  presented  in  Figures  7  and  9. 
Comparison  of  solutions  illustrates  the  excellent  accuracy 
obtained  with  the  Hamming's  Predictor  -Corrector  Method. 


E/E 


pOJ^Q.^ 


FIGURE  4  Correlation  of  Numerical  Integration  Results 
With  Analytical  Solution  For  Centrifigal 
Body  Forces 
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ANALYTICAL  SOLUTION 
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FIGURE  5.  Correlation  of  Numerical  Integration  Results 
With  Analytical  Solution  For  Radial 
Traction  On  Inner  Surface 
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Figure  6.  Correlation  of  Numerical  Integration  Results 
With  Analytical  Solution  For  Radial 
Stress  Traction  On  Outer  Surface 


Table  1  Symbol  Definitions 


01"  SYMBOLS . 


TECi:) 

I  NT 
NN 

Y2LEF'r 

y2r:i:  te 

NHIS 

80:1; 

800 

OTGNL 

DT 

W 

RHO 


initial  starting  ualue  for  independent  variable 

INTEGRATION  STEP-SIZE 


XMAX 

UPPER  LI MI 

N 

NUMBER  OF 

YR  (  I  ) 

N  INITIAL 

FR  ( :i;  ;i 

deriuat  iue: 

EQUATI DNS 

■  ( I V  j ) 

Y  UALUl;::  AT 

■  (  I  ...1 ) 

DERIUAT IUE 

"in  original  SYBl'EM  OF  DIF' FERENTI AL, 

. XH  X  OALUE  FOR  JTH  DIFFERENTIAL  EQUAIION 

/'E  AT  ITH  X  MALUE  FOR  JTH  DIFFERENTIAL 

FOLIATION 

T'  R  0  N  c  A  r  ;i:  0  n  e  r  r  0  f'  0  r  i  t  fi  g  g  R  e;  c  t  0  r  e  0  u  a  1 .1. 0 1'l 

N  0  N  B  E  F';  Li  I--'  I N  T  E  G  R  A  "I  "  1 0  N  8  B  E I W  i::.  F.  N  1)  0  I  i-'  l.l  1  . 

0  p  i"'  E  R  i...  :i:  M  :i:  t  o  n  ’i  h  e  n  d  m  b  e  r  o  i-'  h  a  i...  f  - 1  n  t  f:  i";  v  a  i...  .i.  i  f  i-'  f  J-  o  f:. 

LOWER  LIMIT  ON  UNKNOWN  INITIAL  MALUE 

UPPER  LIMIT  O'N  UNKNOWN  INITIAL  MaLUL  . 

N  U  M  B  E  l''^  0 1"'  H  A I...  I-'  -•  I N  T  E  I'F;  M  A I...  I T I::!  1-5  A 'T 1 0  N  S  B I;:,  i  W  b.  I;:.  N  I-'  R  .1.  N  I  0  U  I 

RADIAL  STRESS  B .  CL  ON  INNER  RADIUS 

RADIAL  STRESS  B .  CL  ON  OUTER  RADIUS 

HAS  THE  UALUE  -•1  v  IF  YKXMAXXSIGO  WHEN 

Y 2 ( X I ) - Y 2 1... if; I" I'  5 0 ri-i i;;; ix w i s e  r i-i i;;:  u a i u f;  i s  i 
'T  E  M  I"’  i;;:  r  a  i"  u  R'  e  d  ;i;  i--  I"  i;;:  fi;  e:  n  c:  f;;  r  e  i...  a  t  :i:  u  f;;  t  o  c:)  j  I’l;  f;;  s  s 

FREE  STATE 

ROTATIONAL  SPEED (RPM) 

Li  i::.  N  S  .1.  I  Y  /  i.  N 


Table  2  Material  Property  Section  of  VARPROP/CFD 


4  ()  0  C 

3  3  ()  ()  C) 
3600  C 
3700 
3800 
3900 
4000 
4100 

4  2  0  0 
4300 
4400  C 
4500  C 
4600  C 
^11^ 


:i:  N  p  u  T  M  A  I'  f;:  FF  ;i:  A  i...  f-'  i-f  o  i-'  f  r  t  i  e;  s 

(  T  i  TANCvE, N  f  I'TI...  v  IX  5  lO  iDIi-'il... )  * 


ET^^Pl,  .4E6 
UTR::::.27 
E:F0:^^2.6E6 
URT"X'F  RYER/ET 
A  T-Pl.  0  1 E-6 

RHO^^^^  <•  I 

f;;  n  d  m  a  t  f;:  f'I!  ;f  a  i...  F'  F';  o  f-'  i;;:  fi;  r  y  d  i;;;  s  c  r  ;f  i-'  t  i  o  n 


Table  3  Material  Property  Section  of  VARPROP/CF  2 


3:1.00  c 

3200.  C 

INPUT  MATER I 

3300  c: 

(  r  t  T  ANGF!NT  :i.  A 

3400  C 

3500 

1::.  1  M  .i.  4  F. 

360  0 

yTR^-.27 

3  7  0  0 

ERM-PI .  SET 

3  &  0  0 

A H  :l.  0  :l.  i:: 

3900 

A  R  M  4  !•  U  F:  - 

0  0  u  L.- 

4 1 0  0  C 

END  MATER  IF 

4200  C 

pROPii-RT  :i:  i:p 

R  j  RAD  ) 


ii...  r’RGF' RRTY  DESCRIPTION 


Table  4  Material  Property  Section  of  VARPROP/NUMD 
(Constant  Properties) 


31 100 'C 

3:1.200  C  INPUT  MATERIAL  PROPERTY  DEF'ENDENCE  ON  RADIAL.  POSITION 
31.300  C 


3 1 4  0  0 

ET^^Pl 

. .  4E6 

3 1 5  0  0 

VTI'O^^ 

-  0  27 

3:1.600 

ER-P 

^ .  6E6 

31700 

U  TRi:: 

*  '  ■  •  0  V 

3  ISO  Qi 

F'l  T  :i 

o.:i.  E-- 

3 .1. 9  0  0 

AR  --"1 

^ ❖ 5E-6 

r'l  1 1 

ERF-^^ 

0  V' 

32:l  00 

ETP- 

(\ 

32200 

A  Tp::; 

0  •> 

32300 

ARI'A^^ 

0  •> 

3  .'•[  4  ( 0* 

ORT- 

ER;};0T 

:.s2500 

OFPri-' 

32600  C  END  HA  TER  I  Al...  PROPERTY  INF'l.lT 

II- 


Table  5  Material  Property  Section  of  VARPROP/NUMD 
(Power  Law  Variation) 


27504  C 

27505  C  ; 

:NPUT  MATEFPCAL  PI'::0PEFPrY 

27506  C 
27600 

ETM^^Pl.  .4E6 

27800 

y-j  R::::  ^  27 

•")  *7  Q  (') 

OTE'P^^^^O. 

28000 

EI-;:M'"4.1 . 6E6 

28006 

OI-PF^^NIRMYOTR/ETM 

:28008 

'v'RTP^^^^O  ■> 

28100 

ATH”"  1 0  V  :!,  E--6 

28300 

A  I'i:  M  4  .1  5 1::!  -•  6 

2831 0 

28312 

ER^^4IRnYXY>f:2 

28314 

AR^-ARii;lXY$2 

28316 

AT-ATM>l<X>Pl:2 

28318 

ERP:^^:2.>  YERMYX 

28319 

ETR^^:^2.KETMYX 

28320 

A  TP -42  .t  ?I<ATM>KX 

2832  1 

ARP^^<P>  YARMI  X 

FND  MATER lAI,.  PROPER  Pv  1 

ON  i-i:AIiIAI,.. 


PGS I T ION 


Table  6  Material  Property  Section  of  VARPROP/NUMD 
(Linear  Variation) 


31:1.00  C 


31200  C 
31300  C 

iNr'i.i  r  MA 

TER  If 

■iL  RRUi-'!:: 

l'^! '  f '  Y  I!l  1::!  [••’  1;;:  N  )J  I:t  N  C  I::!  (!)  N  A  l!  1  A  L. 

POSIT 

31400 

i;:':Ti-T~ 

1 .4E^ 

) 

3:1.500 

OTR::;-' 

A  2  7 

3  :l.  700 

av!  «■  [it.  < 

> 

3  :i.  9  0  0 

UTRI'-' 

-0  « 

3  2  0 0  Q 

A1"M"” 

;!.  0 .  ,1 1 

::  6 

32100 

AF<n- 

4 . 5E” 

••  6 

32200 

ET'-^  ( 

ERM;-I 

FTM  )  YX/  ( 

X  M  A  X  •■■■  X  1  }  ’f  [;•' '  1 "  M  ( 1::!  R  M  - 1::!  1 M  )  :^  X  ] 

/  ( XMA 

32300 

1::!  i’<  C 

ETH-'I 

\RM  )  >y'X/  ( 

X  n  A  X  X  I!  )  1-  FL  F^;  M  -•  ( FF  F  M  •-  F:!  R  M  )  Y  X  FI 

/  ( XMA 

:32400 

AR-  ( ATri"*T 

•iRri )  >I<X/  ( 

X  M  A  X  X  :F  )  -1  A  FL  M  -•  (  A  I"  M  ■■■■  A  R;  M  )  Y  X  FI 

M  A 

32500 

AT-  ( 

A  I'i:  M  A  ’I  M  ,)  :T  X  ( 

X  MAX  ■■■■  X  FIF  )  I-  A  F  M  (  A  IT  M  -•  A  T'  M  )  :=};  X  Fi 

/  ( XMA 

32600 

ERP^::: 

(ETM- 

•TIRM  )  /  ( >; 

M  A  X  X  FIF  ) 

32700 

ETR- 

(ERM- 

Fi:  "F  M  ) ./  ( >■ 

MAX-'X  FIF  ) 

32800 

A 'T  I'" 

( AI-M4' 

••ATM )  /  (> 
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Conclusions 


Development  of  a  numerical  integration  scheme  for 
the  analysis  of  cylindrically  orthotropic  annular  disks 
with  variable  elastic  constants  has  been  accomplished. 

The  integration  scheme  utilizes  Hamming's  Predictor-Corrector 
Method  in  conjunction  with  a  half-interval  search  technique 
which  rapidly  converges  to  the  exact  solution.  Radial 
stress  and  displacement  boundary  conditions  may  be  specified. 
Correlation  of  numerical  integration  results  with  analytical 
and  finite-element  solutions  was  found  to  be  excellent. 

This  analysis  capability  enables  the  determination  of  the 
influence  of  prescribed  material  property  variations, 
temperature  changes,  and  centrifugal  forces  on  the  response 
of  an  annular  disk.  The  disk  may  be  subjected  to  internal 
and  external  pressure  or  displacement  boundary  conditions 
as  well.  In  addition,  interference  fit  can  be  approximated 
by  specifying  displacement  boundary  conditions  on  the 
disk  inner  radius. 
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Appendix  A 

Formulation  of  Governing 
Equations;  Variable  Properties 


For  an  axisymmetrio  body  rotating  at  a  constant 
angular  velocity  w,  the  equilibrium  equation  in  the  radial 
direction  is  given  by; 


2 

pw  r 


0 


(1) 


The  stress  -  strain  relations  in  polar  coordinates  for  a 
cylindrically  orthotropic  material  are  given  as  follows; 


'’e/e 


+  ttg  AT 


(2) 

(3) 


Inverting  (2)  and  (3)  and  solving  in  terms  of  the  stress 
components  yields. 


Or  =  +  Qr0(ee-«0^T) 
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(4) 

(5) 

(6) 


(7) 


The  strain  -  displacement  relations  for  an  axisymmetric 
body  are 


Where  u  is  the  radial  displacement.  To  obtain  the  governing 
equation  in  terms  of  displacement  equations  (4) ,  (5) ,  (8) 

and  (9)  are  substituted  into  the  equilibrium  equation  (1) . 
Simplifying  one  obtains, 
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where 
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and  employing  equations  (6)  and  (7) , 


(19) 


For  the  original  boundary  value  problem,  one  boundary 
condition  (displacement  or  radial  stress  component)  will  be 
prescribed  on  the  inner  and  outer  radii.  Consequently,  the 
unknown  initial  condition  will  be  bounded  and  the  half-interval 
method  will  be  employed  to  iterate  to  the  solution.  Note 
that  Hamming's  Method  always  requires  boundary  conditions  at 
the  inner  radius  in  terms  of  displacement  (see  eq.  (19)). 
Therefore,  if  u(a)  is  prescribed,  then  bounds  on  ^  (a)  are 
established  and  the  solution  is  obtained  in  a  straightforward 
manner.  However,  if  o^{a)  is  prescribed,  bounds  on  00 (a) 
are  expected.  In  this  case,  the  stress  components  are 
transformed  into  the  displacement  and  radial  strain  boundary 
conditions  utilizing  equations  (2)  ,  (.3)  ,  (8)  and  (9)  . 

The  listing  of  the  program  which  performs  the 
numerical  integration  scheme  is  given  below. 
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3100  C  INT  NUMBER  OF  INTEGRATIONS  BETWEEN  OUTPUT 

3200  C  NN  UPPER  LIMIT  ON  THE  NUMBER  OF  HALF-INTERVAL  ITERATIONS 

3300  C  Y2LEFT  LOWER  LIMIT  ON  UNKNOWN  INITIAL  VALUE 

3400  C  Y2RITE  UPPER  LIMIT  ON  UNKNOWN  INITIAL  VALUE 

3500  C  NHIS  NUMBER  OF  HALF -INTERVAL  ITERATIONS  BETWEEN  PRINTOUT 

3600  C  SGI  RADIAL  STRESS  B,C«  ON  INNER-  RADIUS 

3700  C  SGO  RADIAL  STRESS  B^C*  ON  OUTER  RADIUS 

3800  C  SIGNL  HAS  THE  VALUE  -Iv  IF  YKXMAXXSIGO  WHEN 

3900  C  Y2(  XI  )==Y2LEFT  5  OTHERWISE  THE  VALUE  IS  1 

4000  C  DT  TEMPERATURE  DIFFERENCE  RELATIVE  TO  STRESS 

4100  C  FREE  STATE 

4200  C  W  ROTATIONAL  SPEED(RPM) 

4300  C  RHO  DENSITY(#/IN>!<>K3) 

4400  C 

4500  C  SUBROUTINES  RUNGE  AND  HAMMING  ARE  BASED  ON  ALGORITHM'S 

4600  C  FOUND  IN  "APPLIED  NUMERICAL  METHODS"  BY  CARNAHAN ^  LUTHER 

4700  C  AND  WILKES  (PAGES  367-402) 

4800  C 

4900  C  MAIN  PROGRAM 

5000  C 

5100  INTEGER  COUNT v RUNGE r HAMING 

5200  LOGICAL  PRED 

5300  DIMENSION  TE(IO) c YR ( 10 ) y FR ( 10 ) y Y ( 4 y 10 ) » F< 3 r 10 ) y PHI (10) ySAVEY(lO) 

5400  lyYPRED(lO) 

5500  COMMON  DT  y  W  y  RHO  y  X I y  XMAX 

5600  C 


READ  INPUT  DATA 


atiw  i, 
5900 
6000 
6100 
6200 
6300 
6400 
6500 
6600 
6700 
6800 
6900 
7000  C 
7100  C 
7200  C 
7300 
7400 
7500 
7600 
7700 
7800  C 
7900  C 
8000  C 
8100  C 
8200  C 
8300 
8400 
8500 
8600 
8700 
8800 
8900  C 
9000  C 
9100  C 
9200 
9300 
9400 
9500 
9600 
9700 
9800 
9900  C 
10000  C 
10100  C 
10200  C 
10300  C 
10400 
10500 
10600 
10700  C 
10800  C 


N=2 

WRITCi;<6rlO) 

READ  ( 5 / )  XI  y  H  V  XMAX » INT y  DT  y  W y  RHO 

WR;rrE(6y20) 

READ(5y/)YRlyTlCl 

WRITE(6y25) 

READ ( 5  y / ) Y2LEFT  y  Y2R I TE 
WRlTE(6y26) 

READ(5y/)YR3yTlC3 
WRITE (6y 27) 

READ ( 5  y / ) NN  y  NH I S  y  SI GNL 


F’RINT  HEADING  AND  INF‘'UT  DATA 


WRITE ( 6  y  30 ) Hy  XI y  XMAXy  N  y INT  y  YRl y  TICl y  Y2LEFT  y  Y2RITE  y  NN  y  NHIS  y  SIGNL 
1  yDTyRHOyW 
WRITE(6y40) 

W===W*2.!{c3*  1415927/60. 

RH0:==RH0/32.2/12. 

INITIALIZE  STEP  COUNTER 

SET  FIRST  ROW  OF  Y  MATRIX  EQUAL  TO  INITIAL  VALUES 
INITIALIZE  TRUNCATION  ERRORS  TO  ZERO 

DO  90  II"lyNN 
X:==XI 

YR(1)-:YR1 

Y2ZER0'-=  ( Y2LEFT  f  Y2R I  TE )  /2  ♦ 

YR(2):=i:Y2ZER0 

WRITE ( 6  y  35 ) 1 1 y  Y2LEFT  y  Y2RITE  y  Y2ZER0 
CONVERT  STRESS  B.C.  TO  DISPLACEMENT  B.C. 

IF ( TICl. EQ. 2) GO  TO  55 

CALL  F'Ia-OP  (  X  y  YR  ( 1 )  y  YR  ( 2 )  y  FF  y  P  y  Q  y  YR  ( 2  )  y  EPT  y  YR  ( 1  )  y  1) 

55  COUNT==:0. 

M^=0. 

DO  60  J=lyN 
TE(J)^==0. 

60  Y(4y  J)"-=YR(  J) 


CALL  FOURTH -ORDER  RUNGE-KUTTA  FUNCTION  TO  INTEGRATE 
OVER  FIRST  THREE  STEPS 

SUBROUTINE  DERIV  CALUCULATES  DERIVATIVES 

65  I F ( RUNGE ( M  y  N  y  YR  y  FR  y  X  y  H  y  PH I y  SA VEY ) . NE . 1 ) GOTO  70 
CALL  DERIV ( N  y 1 y ISUB  y  YR  y  Y  y  FR  y  F  y  X ) 

GOTO  65 


PUT  APPROPRIATE  INITIAL  VALUES  IN  Y  AND  F.  MATRICES 


•••» 


/•%/’%!  I  >  I  : _ /•\  /•\t  t\t’f  I 


XV'fW  /U  UUUIV  I  “UUUIN  I  t  .1. 

11.000  i:SUB==4“C0UNT 

11100  HO  75  J==li-N 

11200  75  Y(3;SUB.J)==YF«(J) 

11300  CALL  BERI0(iM>2!-;i;SUB»YRfYi-FR»FrX) 

11400  C 

11*500  c  pr;i;nt  solut;i;ons  after  int  steps 

11600  C 

11700  80  :IF(C0UNT/.INT*;i;NT.NE*C0UNT)G0  to  85 

1 1  BOO  .i; F  <  .  NOT  *  ( 1 1 /NH 1  S>KNH I S ,  EQ »  X I .  OR .  1 .1  ♦  EO  *  NN )  )  GO  TO  85 

11900  lF(C0UNT.LE.3)EPT=Y<:i;SUBvl.  )/X 

12000  IF(  COUNT ♦LE*3)CALL  PROP  ( X  >  F(  iSUBv  1 )  s-EPTrFFyPj  Oi-SIGRr  SIGT  y  Uf  3  > 

12100  IF(C0UNT.GT«3)EPT:=Y<lyl)/X 

12200  IF ( COUNT ♦ GT ♦ 3 ) CALL  PROP ( X  >  F ( 1 y 1 ) y  EPT  y  FF y  P  y  0  y  SIGR  y  SIGT  y  U  y  3 ) 

12300  IF(C0UNT«LE*3 ) WRITE (6 y 50 ) X  y  Y ( ISUBy 1 ) y F( ISUBy 1 ) y EPT y SIGR  y SIGT 

12400  IF ( COUNT ♦ GT ♦ 3 ) WRITE ( 6  y  50 ) X  y  Y ( 1 y 1 ) y  F ( 1 y 1 ) y  EPT  y  S IGR  y  SIGT 

1 2500  C 

12600  C  IF  X  >  XMAX  TERMINATE  INTEGRATION 

12700  C 

12800  85  CONTINUE 

12900  IF<X*GT*XMAX-H/2)G0T0  100 

13000  C 

13100  C  CALL  RUNGE  OR  HAMMING  TO  INTEGRATE  NEXT  STEP 

1 3200  C 

13300  IF ( COUNT »LT. 3) GOTO  65 

13400  PRED-~«TRUE* 

13500  105  MM~I  IAMING  ( N  y  Y  y  F  y  X  y  H  y  TE  y  PRED  y  YPRED  ) 

1 3600  CALL  DER I M ( N  y  3  y I SUB  y  YR  y  Y  y  FR  y  F  y  X ) 

13700  I F  <  MM ♦ EO ♦ 1 ) GOTO  1 05 

13800  C 

13900  C  INCREMENT  STEP  COUNTER  AND  CONTINUE  INTEGRATION 

14000  COUNT:=COUNTH 

14100  GO  TO  80 

14200  C 

14300  C  COMPARE  SOLUTION  TO  KNOWN  OUTER  B,C* 

14400  C 

14500  100  IF<TIC3.EQ.2)G0  TO  102 

14600  C 

14700  C  CONCERT  DISP*  AND  STRAIN  TO  STRESS 

14800  C 

14900  EPT==“Y(lyl)/X 

15000  EPR=“H"' <  1  y  1 ) 

15100  CALL  PR0P(XyEPRyEPTyFFyPy0ySIGRySIGTyUy3) 

15200  IF(  (SIGR-YR3)>KSIGNL*6T,0>G0  TO  106 

15300  Y2RITE===Y2ZERG 

15400  GO  TO  90 

15500  106  Y2LEFT====Y2ZER0 

15600  GO  TO  90 

15700  102  IF(  (Y(lyl)“YR3)>KSIGNL  ♦GT.  0)G0  TO  110 

15800  Y2RITE==Y2ZER0 

15900  GO  TO  90 


16000 

16100 

16200 

16300 

16400 

16500 

16600 

16700 

16800 

16900 

17000 

17100 

17200 

17300 

17400 

17500 

17600 

17700 

17800 

17900 

18000 

18100 

18200 

18300 

18400 


110  Y2LEFT-Y2ZERQ 
90  CONTINUE 

FORMAT  STATEMENTS 

10  FORMAT  ( IX  y  '  ENTER  XI  y  H  y  XMAX  y  INT  y  DT  y  W  ( RPM )  y  RHO  (  LB/IN>K>K3  )  '  ) 

2 0  I"  0 1^ M A T  ( 1 X y  ' E N T E R  K N 0 Ul N  I »  y  I’ Y P E (  S T R E S S -=  1  y  D I S P L  » 2)  '  ) 

25  FORMAT  < 1 X  y ' ENTER  Y2LEFT  y  Y2R I TE ' ) 

2 6  I"' 0 R M A  T ( 1 X y  ' E N I' E R  S E C 0 N B  B » C ♦  y  T Y P E ( S T R E S S  ==•  1  y  B I S P L  * ••=  2)  '  ) 

2  7  I"  0 1^  M  A  T  ( 1 X  y  '  E  N  T  E  R  N  N  y  N  M I S  y  S I G  N I... '  ) 

30  FORMAT ( ///IX y  ' H  '  y  E15 ♦  3/lX y  ' XI  "" '  y  F15 ♦  4/ 

J  IXy'XMAX  yF15*4/lXy 'N  == '  y  1 15/lX y  MNT  =-'yI15/ 

1  IX  y  '  YRl  - '  y  F15 , 4/lX  y  '  TICl  == '  y  FIS  ♦  4/lX  y  '  Y2LEFT== '  y  F15  ♦  4/ 
1  IXy  '  Y2RITE“'' yF15*4/lXy 'NN  y  IlS'/lX^NHIS  ='yI15/ 


y FI 5*4/1  Xy  "RHO 


yF15*4/lXy 


1  IXy 'SIGNL.  ==='I15/lXy 'BT  =  '  y  F15  *  4/lX  y  "  RHO  =  '  y  F- 15  ♦  4/lX  y 
1  "W  ::="F15*4//) 

35  FORMAT  ( 1 X  y  "  I  TER  AT  I  ON "  y  1 5  y  5X  y  '  LOWER=== '  F 1 5  ♦  4  y  5X  y  '  UPPER=  '  y  F 1 5 . 4 
1  5X  y  "  UNKNOWN  I  *  C  *  '  F 1 5  ♦  4 ) 

40  FORMAT ( 7X  y ' X ' y  30X  y ' Y ' y /30X  y ' U " y 18X  y ' EPR " y 16X  y ' EPT " y 16X 
1  y "SIGR" yl6Xy "SIGT"/) 

50  FORMAT ( 1 X  y  F 1 0  *  5  y  5X  y  5 ( 5X  y  E 1 5 ♦ 7 ) ) 

CALL  EXIT 
ENB 

C  )K  )f<  :<<>!<  >K  *  )K  >{«>!;  JK  $  >l<  >1<  >!<  >K  *  JK*  !f{  ^:  >!«  >K  >K  :<<  $  sK*  >K  >K  >|c  ^  )K ){!)!«:{{ JK  *  !F;  5{{  )1«  )K  )1(  5l<  ^  >K  )!<*****)!<>!< 

F  U  N  C  I’  1 0  N  R  U  N  G  E  ( M  y  N  y  Y  y  F  y  X  y  I- 1  y  P  H I  y  S  A  V  E  Y ) 


18500  C 

18600  INTEGER  RUNGE 

18700  B I MENS I ON  PH I ( N ) y  S A VE Y ( N ) y  Y  <  N ) y  F ( N ) 

18800  M==MH 

1 8900  GO  TO  ( 1 y  2  y  3  y  4  y  5 ) y  M 

19000  C 

19100  1  RUNGE" 1 

19200  RETURN 


19300 

C 

19400 

2 

BO  22  J==-lyN 

19500 

SAVEY( J)"Y( J) 

19600 

PHI  ( J)^==F(  J) 

19700 

Y  (  J  )  "SAVEY  ( J )  i-  *  5>KH*F  (  J 

19800 

X”X{“*55i?H 

19900 

RUNGE==1 

20000 

RETURN 

20100 

C 

20200 

3 

BO  33  J:-=lyN 

20300 

PHK J)"PHI<J)i2*^F( J) 

20400 

33 

Y  ( J )  "-^TSAOEY  ( J )  ■{•  *  5)i<H>!<F  <  J ) 

20500 

RUNGE~Ll. 

20600 

RETURN 

20700 

C 

20800 

4 

BO  44  J=:=vlyN 

20900 

PHK  J)=^PHI(J)i2*>KF(  J) 

21000 

44 

Y  ( J )  ==SAMEY  ( J )  }  H>i<F  (  J ) 

2i:i.oy 

21200 

2:1.300 

21400 

21500 

21600 

21700 

2.1‘aoo 

21900 

22000 

22100 

22200 

22300 

22400 

22500 

22600 

22700 

22800 

22900 

23000 

23100 

23200 

23300 

23400 

23500 

23600 

23700 

23800 

23900 

24000 

24100 

24200 

24300 

24400 

24500 

24600 

24700 

24800 

24900 

25000 

25100 

25200 

25300 

25400 

25500 

25600 

25700 

25800 

25900 

26000 

26100 


X“Xi»5)KH 

Fi:UNGE="l 

RETURN 

5"  DO  55  J~lyN 

55  Y  ( J )  ^SAOE  Y  <  J )  f  <  PI  1 1  ( J )  f  E  ( J )  )  >KH/6 

M=0. 

RUNGE=--0 

RETURN 

END 


SUBROUTINE  DERIO ( N r FCOUNT v ISUB » YR y Y y FR y F y X ) 

DIMENSION  YR<N) y Y(4yN) yFR(N) yF(3yN) 

CALL  PROP ( X  y  S I GR  y  S I GT  y  FF  y  P  y  0  y  EPR  y  EPT  y  U  y  2 ) 

GO  T0(ly2y3)y  FCOUNT 

1  FR<1)=YR(2) 

FR  <  2 )  =FF"-P*YR  ( 2 ) -QtYR  <  1 ) 

RETURN 

2  F(ISUByl)==YR(2) 

F  (  I  SUB  y  2  )  =FF-p>K  YR  ( 2 )  -a>KYR  ( 1 ) 

RETURN 

3  F(lyl)==Y(ly2) 

F  ( 1  y  2 )  =FF-  P)I<Y  ( 1  y  2 )  "“Q^Y  ( 1  y  1 ) 

RETURN 

END 


FUNCT ION  MAM ING ( N  y  Y  y  F  y  X  y  H  y  TE  y  PRED  y  YPRED ) 

INTEGER  NAMING 
LOGICAL  PRED 

D I MENS I ON  YPRED ( N ) y  TE ( N ) y  Y ( 4  y  N ) y  F ( 3  y  N ) 

IS  CALL  FOR  PREDICTION  OR  CORRECTOR  SECTION 
IF < y NOT. PRED)  GOTO  4 

PREDICTOR  SECTION  OF  NAMING 
COMPUTE  PREDICTED  VALUES  AT  NEXT  POINT 
DO  1  J=lyN 

YPRED  (  J )  ~Y  ( 4  y  J )  •1-4 .  *N>K  ( 2  ♦  >!«F  ( 1  y  J )  ~F  <  2  y  J )  12  *  ( 3  y  J  )  )  /3 . 

UPDATE  THE  Y  AND  F  TABLES 
DO  2  J-lyN 
DO  2  K5=-=ly3 
K=5"K5 

Y(KyJ)===Y(K"-lyJ) 

IF ( K . LT . 4 ) F ( K  y  J ) ™F ( K-1 y  J ) 

CONTINUE 

MODIFY  PREDICTED  Y(J)  VALUES  USING  THE  TRUNCATION  ERROR 


2620U  i; 
26300 
26400 
26500 
26600  C 
26700  C 
26800 
26900 
27000 
27100  C 
27200  C 
27300  C 
27400  C 
27500 
27600 
27700 
27800 
27900  C 
28000  C 
28100 
28200 
28300 
28400 
28500  CX 
28600  C 
28700 
28800  C 
28900  C 
29000  C 
29100  C 
29200  C 
29300  C 
29400  C 
29500  C 
29600  C 
29700  C 
29800  C 
29900  C 
30000  C 
30100  C 
30200  C 
30300  C 
30400  C 
30500  C 
30600  C 
30700  C 
30800  C 
30900 
31000 
31100  C 
31200  C 


hyilflAlki-S  l-KUn  nil;:.  I 'KhV  l-UU;;)  5  I  tl  •  }•  .lNCKI;:.ril;:.N  I  X  VHLUt. 

no  3  J-li-N 

3  Y  <  1  y  J )  =  YF'REH  ( J )  H 1 2 » >i<TE  ( J )  /9 , 

X=X-fH 

SET  PREn  AND  REQUEST  UPDATED  DERIVATIVE  VALUES 
PRED==  ♦  FALSE . 

HAMING=1 

RETURN 


4 


CORRECTOR  SECTION  OF  HAMING 

COMPUTE  CORRECTED  AND  IMPROVED  VALUES  OF  THE  Y(J)  AND 
SAVE  TRUNCATION  ERROR  ESTIMATES  FOR  THE  CURRENT  STEP 


DO  5  J=lyN 

Y(  J  y  J)  =  (9.>i<Y<2y  J)“  Y(4y  J)i3.5KH>f«<F(lyJ)+2 

TF  ( J )  ==9 ,  ( Y  <  1  y  J )  -YPRED  (  J  )  )  / 1 2 1  ♦ 


)i{F(2y  J)“F(3y  J)  )  )/8. 


Ydy  J)==Y(lyJ)“TE(J) 


SET  PRED  AND  RETURN  WITH  SOLUTIONS  FOR  CURRENT 
PRED====  >  TRUE  ♦ 

HAMING==2 

RETURN 


SUBROUT I NE  PROP ( X  y  A 1 y  A2  y  FF  y  P  y  Q  y  A3  y  A4  y  A5  y  OPT 1 ) 


INPUT  MATERIAL  PROPERTY  DEPENDENCE  ON  RADIAL  POSITION  X 
( T i  TANGENT I AL  y  R  t RAD I AL ) i 


ET  -TANGENTIAL  MODULI 
FTP  -DERIVATIVE  OF  ET 
V  T  R  -  P  0 1 S  S  0  N  R  A  T 1 0 
VTRP  -DERIVATIVE  OF  VTR 
ER  -RADIAL  MODULI 
ERP  -DERIVATIVE  OF  ER 
VTR  -POISSON  RATIO 
VTRP  -DERIVATIVE  OF  VTR 

AT  -TANGENTIAL  COEFFICIENT  OF  THERMAL  EXPANSION 
ATP  -DERIVATIVE  OF  AT 

AR  -RADIAL  COEFFICIENT  OF  THERMAL  EXPANSION  _ 

DT  -TEMPERATURE  Cl lANGE (POSITIVE  VALUE  COKKOSPONDS  10  AN  INLKEAaL 
RELATIVE  TO  THE  STRESS  FREE  STATE) 

EPR  -RADIAL  STRAIN  COMPONENT 

EPT  -TANGENTIAL  STRAIN  COMPONENT 
U  -RADIAL  DISPLACEMENT 


REAL  K 

COMMON  DTyWyRHOyXI yXMAX 

INPUT  MATERIAL  PROPERTY  DEPENDENCE  ON  RADIAL  POSITION 


31.  >500  C 
31400 
31500 
31700 
3iyoo 
32000 
32100 
32200 
32300 
32400 
32500 
32600 
32700 
32800 
32900 
32950 
32960 
33000  C 
33100 
33200 
33300 
33400 
33500 
33600 
33700 
33800 
33900 
34000 
34100 
34200 
34300 
34400 
34500 
34600 
34700 
34800 
34900 
35000 
35100 
35200 
# 


E:TM==1  *4116 
VTR=-“- » 27 


I'  ••••  0  r  R  E  R  /  E  T  >i<  2  E  T  P 
RTY  INPUT 


ERM^=^2*6E6 
VRTP^^O. 

ATM-IO. lE-6 

E  T  ( E  R  M  ••••  E  T  M )  X  /  ( X  M  A  X  "  X I )  f  E  T  M  -  ( E  R  M  •  •  E  I’  M )  >f<  X I  /  ( X  M  A  X  -  X  .1. ) 
ER-"  ( ETM  "  ERM )  >KX/  ( XMAX  -XI )  -f  ERM  -  ( ETM-ERM )  )ifXI/  ( XMAX-XI ) 
AR""  ( ATH-  ARM )  tY./ <  XMAX-X I )  iARM-  ( AIM-ARM )  >I<X  1  /  ( "  wi  \ 
AT-  ( ARM -ATM )  >KX/ ( XMAX-XI )  +ATM-  ( ARM-ATM ) >i<XI/  ( XMAX-XI ) 

E  R  P  ===  <  E  r  M  -  E  R  M )  /  <  X  M  A  X  -  X I ) 

E  I’  P  ==  ( E  R  M  •••■  E  T  M )  /  ( X  M  A  X  -  XI) 

ATP” ( ARM-ATM ) / ( XMAX-XI ) 

ARP===  <  ATM-ARM ) / <  XMAX-XI ) 

ORT^=ER>!<VTR/ET 

V  R  r  p  -  0  T  R  *  E  R  p  / 1;-  T  •••■  0  r  r  t  e  r  /  e  t  C:.  i  p 

MATERIAL  PROPERTY  INPUT 
VI:m:=1  .-VRT:>1«VTR 
QRR""0::;R./VD 
QTT::^=ET/Vi:i 

QRT:^=ER>f:OTR/OD  . . 

(!■!  F?  F<  F-’  “=  F:  F^  F-'  /  Ei!  F^  }■  0  F<  I' 1 0  T  F?  F-'  /  0  D  T  0 1 FH  ^  V  F\  F  P  /  L  D 
Q  F^  T  P  ••=  0  r  Fi;  E:!  R  F‘  /  E:!  Fi:  •}■  ( 0  T 1'^  P  i'  'v'  T  Fi:  5j<  !K  2  $  0  F^  T  F  ' )  /  V  D 
K::=SaRT(ET/ER) 
if<optij:q»  1  )G0  to  i 
I F"  ( 0  P  T 1  *  E  0 . 3 )  0  0  7  0  3 

I"  p  -  F?  FF  0  5f'  W  >l«  $  2  >!<  X  /  Q  I-?  Fi!  f  ( <  (7  R  R  P  f '  <  1  0 1 R )  /  X )  ¥  A  F< 

( Q  R  T  F-'  •}■  K  ¥  ¥  2  ¥  ( 0  R  T  - 1  ♦  )  /  X )  ¥  A  T  f  A  R  P  F-  V  7'  R  ¥  A  7'  P )  ¥  H  7 
p:=l./X}QRRP 
0:---=aRTP/X-(K/X)>!<>K2 

if(optij;;;q.2)go  to  2 
A  3  ==  A 1  /  E  R  -  V  r  Fi:  ¥02/  E I  f  A  R  ¥  D  7' 

A  4  ;=  -  V  Fi:  7'  ¥  A 1  /  E  R  •}•  A  2  /  E  7'  f  A  T  ¥  D  T 
A5™X¥A4 
00  70  2 

A3==T;!Fi:Fi:¥A  1  f  QFi:T¥A2-i;iT¥  ( QFi:Fi:¥ARiQFi:T¥A7  ) 

A4-0Fi:T¥Alf  QTT¥A2-DT¥  ( OFi:T¥AFi:}OTT¥AT ) 

RETURN 
END 


^  A  Egd+kv^j) 


‘i-Wer* 


‘^"^re^0r^ 


(l+3v^.)k^pw^r^  ATEg 


If  radial  stress  boundary  conditions  are  prescribed 
on  the  inner  (a)  and  outer  radius  (b) , 


(a)  =  p 


cT^(b)  =  q 


The  unknown  constants  are , 


-  Q) 


bk-l[d2k_i]  Q^^(k+Vg^) 


(-P  +  d^~^Q)d^^^ 
b-k-ljd^’^-l]  Qj.j.(V0r-k) 


where  P  =  p 


<3+'>er>DioV  ATEg(a^-ag) 
Q  =q+  — ——2  -  2 


d  =  a/b,  k^  =  Eg/E^ 


Employing  equation  (26)  ,  we  obtain  the  following  relations 


for  the  stress  components : 


(28) 


(29) 


A  program  listing  for  the  analytical  solution  derived  for 
radial  stress  boundary  conditions  is  given  below.  Stress 
components  are  determined  directly  from  equations  (28)  and 
(29)  and  strain  components  and  radial  displacement  are  calculated 
from  equations  (2),  (3)  and  (9),  respectively. 


#F.tLE 

1000 

1100 

1200 

1300 

moo 

1500 

1600 

1700 

1800 

1900 

2000 

2100 

2300 

2500 

2600 

2800 

2900 

3000 

3100 

3200 

3300 

3350 

3360 

3400 

3500 

3600 

3700 

3800 

4000 

4200 

4250 

4300 

4500 

4550 

4600 

4700 

4800 

4900 

5000 

5100 

5200 

5300 

5400 

5500 

5600 

5700 

5800 

5900 


( 1345)0ARPRaP/CFri  ON  PACK 


$RESET  FREE 

FILE  6  (KIN  D  ==  R  E  M  0  T  E  y  M  A  X  R  E  C  SIZE  2  2 ) 

REAL.  K 

c  LIST  OF  symbols: 

C 

C  XI  INITIAL  STARTING  VALUE  FOR  INDEPENDENT 


C 

C 

C 

C 


DX  X  INCREMENT 

XMAX  UPPER  LIMIT  OF  INTEGRATION 
SGI  RADIAL  STRESS  B.C*  ON  INNER  RADIUS 

SGO  RADIAL  STRESS  B*C,  ON  OUTER  RADIUS 


VARIABLE 


C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 


ET  -TANGENTIAL  MODULI 

VTR  -POISSON  RATIO 

ER  -RADIAL  MODULI 

AT  -TANGENTIAL  COEFFICIENT  OF  THERMAL  EXPANSION 
AR  -RADIAL  COEFFICIENT  OF  THERMAL  EXPANSION 

D T  - T E M P E R A  T U I’v' E  C H A N G E  ( I-'' 0 S 1 1'  I V E  V A I... U E  C 0 R R 0 S l■•' 0 N D S  T 0  AN  INCREASE 

RELATIVE  TO  THE  STRESS  FREE  STATE) 

EPR  -RADIAL  STRAIN  COMPONENT 
EPT  -TANGENTIAL  STRAIN  COMPONENT 
U  -RADIAL  DISPLACEMENT 

W  -  R  0  T  A  T 1 0  N  A I...  S  P  E  E  D  ( R  P  M ) 

RHO  -DENSITY  (#/IN>i<>!<3) 


C 

C  INPUT  MATERIAL  PROPERTIES 
C  (T:  tangential y  R: RADIAL): 

C 

ET===2*6E6 

VTR=-=«5()1 

ER=^imE6 

vrt==vtr;ker/et 

AT:^=4.5E-6 

AR==10«lE-6 

RHO==.l 

C 

C  END  MATERIAL  PROPERTY  DESCRIPTION 
C 

C  READ  INPUT  DATA 

C 

WRITE(6y:l.0) 

READ  (  5  y ./  )  X I  y  XMAX  y  DX  y  DT  y  W 
WRITE (6 y 20) 

READ(5y/)SGI ySGO 
C 

C  PRINT  HEADING  AND  INPUT  DATA 

C 

W  R'  I T  E  (  6  y  3  0  )  X I  y  X  M  A  X  y  S  G  0  y  S  G I  y  E  1"  y  E  R  y  A  T  y  A  R  y  V  I’  R  y  D  I'  y  W  y  R 1 1 0 
WRITE(6y40) 


y.  fs 


W==W*2>K3»  1415927/60* 

RHa=RH0/32*2/12. 

D=X.i:/XMAX 

K=:=SaRT(ET/ER) 

liEN~D)!<)l{(2>lcK)-l. 

PI  ==IiT>fc ( AR“  AT ) / ( 1 .  •-K>!<^2 ) 

P2===  ( 3iMTR )  >I<RI  IO>l<W^'(>}c2/  ( 9 .  ~K)K.lf2 ) 

P = S  G I  f  P  2  )K  X I  Xc  1 2  -  P 1  )K  K  X:  X:  2  $  E  R 

a=-SGC)+P2X<XMAXXcX{2"-PlXcKX<X«2X<ER 

X=XI 

R==X/XMAX 

S;iGR=<PX{i:iX<X<(Kfl)-Q)X<RXcX<(K-l  )/DEN~(P"DX:X<(K-"l  )Xta)X:DX'Xc(K+l ) 
/DENXcRX<;K  ( -K- 1 )  “P2XcXX:X(2-f  PI  X:KX<X<2X:ER 
SIGT"KXcXc2X<  ( I  fKXcORT )  /  ( K  f  VTR )  t  (  PXcDjKXc  ( K+1 )  -  Q )  /DENXcRXcX'  ( K-1 ) 
-KX<X<2Xc  ( 1  -KXcORT )  /  ( OTR-K )  X<  ( P-riX<X(  ( K- 1 )  XcQ )  /  DENXiUXtXc  ( Ki  1 ) 

X«RX<X<  ( -K- 1 )  -  ( 1 +3  *  X«MRT )  X!RHOXc  ( WX<ia'X )  X<X<2/  ( 9-Ktt2 )  +P 1  X<ET 
EPR==S  I  GR/ER-OTRXtS  1  GT/ETIARXcDT 
EPT:=-0TRX<S  I  GR/ETiS  I GT/ET+ ATX<DT 
U==X:XfEPT 

W  R I T  E  ( 6  y  5  0 )  X S 1 G  R  y  S I G  T » E  P  R  i-  E  P  T  r  U 
IF(X*GE*XMAX-DX)GO  TO  2 
X=XfDX 
GO  TO  80 

FORMAT  STATEMENTS 

FORMAT < IX y " ENTER  XI y XMAX y DX y DT y W " ) 

FORMAT ( 1 X  y ' ENTER  SGI y  SGO ' ) 

FORMAT(///lXy  "XI  :== "  y  E15  *  3/lX  y  "  XMAX  ='fE15*A/ 

IXy'SGO  ===' yE15.4/lXy  "SGI  === "  y  El  5  *  4/1 X  y  "  ET  *:"yE15.4/ 

lXy"ER  ••="  yE15.4/lXy  "AT  ="  y  E15  *  4/lX  "  AR  ~"yE15*4/ 

lXy"OTR  ="  yE15*3/lXy  "CT  ==  "  y  E 15  ♦  4/lX  y  "  W  ==="yE15*4/ 

lXy"RHO  ="E15*4//) 

FORMAT(  17Xy  "X"  y  14Xy  "SIGR"  y  16Xy  "SIGT"  y  16Xy  "EF'R"  y  17X 
y "EPT" yieXy "U"/) 

FORMATS  10XyF10.4y5(5XyE15*7) ) 

RETURN 

END 


Appendix  C 


Analytical  Solution:  Power  Law  Variation  of  Properties 


For  the  power  law  variation  of  modulii,  the  stress 
function  approach  [4]  is  utilized  to  obtain  a  closed  form 
solution.  Employing  the  strain-displacement  relations  in  (8) 
and  (9)  and  eliminating  u  in  (2)  and  (3) ,  we  obtain 


^^Gr  ^r  r 

(  -  +  (I  “e)  + 

r  0 


^  -  IT-  ^6 

r  6 


where  prime  denotes  derivates  with  respect  to  radial  position.  , 
For  an  axisymmetric  body,  the  relations  between  the  stress 
function  and  stress  components  which  identically  satisfy 
Eq.  (1) ,  (in  the  absence  of  body  forces)  are 


Oj.  =  g/r,  Og  =  g' 


Substitution  of  (26)  into  (25)  yields  the  governing  differential 
equation. 


I  E.'  ^^fiyEfl' 

g"  +  {-  -  _l_)a'  +  (  ^  ■ 

^  r  E  ^g  rE 

tg  r£.g 


=  -  AT  [ 


E  (a  -a  ) 

e  e  r' 


+  ] 

e  0 


Assuming  the  exponential  form  for  the  modulii  equation  (32) 
may  be  ordered 


E  = 
r 


=  ^rm  ^  '  ^0=  ^m 


=  constant,  =  Vg^  ^ 


r^g"  +  r(l-m)g'  +  (mVg^-k^)g 


^rr,  T?  f  ^*^6  °^r^  ,  ,  , 

-  AT  Eg^  r  [  - ^ -  +  ae'] 


The  general  solution  of  equation  (29)  can  be  expressed  as. 


N  N 

g  =  Ar  +  Br  +  g. 


where  N, 


m+  m^  +  4(k^-inUg^) 


and  is  the  particular  solution. 

It 

The  particular  solution  is  determined  by  a  variation 
of  parameter's  approach, 

Ni  N2 

gp  =  v^r  +  V2r 


where 


-r  C. 
^l'  =  ■(N2-N1) ' 


1-N« 

^2'  (N2-N^) 


and  is  defined  as  the  right-hand  side  of  equation  (34) 


^  m+2  ,  ,  ,, 

■^3  =  f —  +  “e  ^ 


The  stress  components  defined  in  equation  (31)  are; 


Nf-l  N2-I  N,-l  Np-l 

=  Ar  +  Br  +  r  +  r  V2 


N,-l  N2-I  N,-l  N2-I 

Oq  =  r  +  BN2  r  +  N^^r  +  N2r 


and  A  and  B  are  unknown  constants  to  be  determined  from  the 
boundary  conditions.  The  components  of  displacement  and 
strain  are  calculated  from  equations  (9),  (2)  and  (3), 

respectively. 

If  radial  stress  boundary  conditions  are  prescribed 
in  equatiop  (25),  the  unknown  constants  are  given  as  follows; 


N,-l  N,-N, 
b  (d  -1) 


N,-l  l-N- 
-  (P-d  Q)d 
Np-l  N,-N2 
b  ^  (d  ^-1) 


where 


N,-l  N2-I 

P  =  p  -  a  V^(a)  -  a  ^  V2(a) 


N,-l  N^-l 

Q  =  q  -  b  Vj^(b)  -  b  ^2  (b) 


and  the  (i=l,2)  are  obtained  through  integration  of 
equation  (37)  and  evaluation  at  the  boundaries. 

For  completeness,  assume  the  coefficients  of  thermal 
expansion  also  vary  as  a  power  law,  i.e.. 


=  «em  ^ 


a  =  a  r 
r  rm 


Equation  (38)  reduces  to 


C3  -  c^r 


2m"l 


where 


C4  = 


and  integrating  equation  (37)  yields 


-C4r 


2  m-Nj^+1 


(N2-N1)  (2m-N^+l) 


2m  -N,+l 

V,  =  _ 

2  (N2-Nj_)  (2m-N2+l) 


Employing  equations  (41),  (42),  (43)  and  (44)  we 
obtain  the  following  expression  for  the  stress  components: 


Note  that  for  m  =  0  (Nj^=k,  N2=-k)  ,  equations  (51)  -  (  54) 
reduce  to  the  solution  given  in  Appendix  B  for  constant 
material  properties.  The  program  listing  for  this  solution 
is  given  next. 


#FILE  <1345)VARPR0P/CF2  ON  PACK 
1000  PRESET  FREE 


1100  FILE 
1200 
1300  C 
1400  C 
1500  C 
1600  C 
1700  C 
1800  C 
1900  C 
2000  C 


6  ( K I  ND==REMOTE  y  MAXRECS I  ZE~22 ) 

REAL  KyNlyN2yM 

LIST  OF  symbols: 

XI  INITIAL  STARTING  VALUE  FOR  INBEPENBENT  VARIABLE 

BX  X  INCREMENT 

XMAX  UPPER  LIMIT  OF  INTEGRATION 
SGI  RABIAL  STRESS  B«C.  ON  INNER  RABIUS 

SGO  RABIAL  STRESS  B.C.  ON  OUTER  RABIUS 


2100  C 
2200  C 
2300  C 
2400  C 
2500  C 
2600  C 
2700  C 
2800  C 
2900  C 
3000  C 
3100  C 
3200  C 
3300  C 
3400  C 
3500 
3600 
3700 
3800 
3900 


ETM  -TANGENTIAL  MOBULI 
VTR  -POISSON  RATIO 
ERM  -RABIAL  MOBULI 

ATM  -TANGENTIAL  COEFFICIENT  OF  THERMAL  EXPANSION 
ARM  -RABIAL  COEFFICIENT  OF  THERMAL  EXPANSION 

BT  “TEMPERATURE  Cl lANGE (POSITIVE  VALUE  CORROSPONBS  TO  AN  INCREASE 
RELATIVE  TO  THE  STRESS  FREE  STATE) 


EPR 

-RABIAL 

STRAIN 

COMPONENT 

EPT 

-TANGEN' 

riAL  STI 

a'AIn  component 

U 

“RABIAL 

BISPLACEMENT 

INPUT 

MATERIAL 

PROPER 

riES 

(T:  tangential y  R: RABIAL): 


ETM==1*4E6 
VTR==»27 
ERM==2 . 6E6 
ATM=10*1E“6 
ARM=4 . 5E-6 


4000  C 
4100  C 
4200  C 
4300  C 
4400  C 
4500 
4600 
4700 
4800 
4900  C 
5000  C 
5100  C 
5200 
5300 


ENB  MATERIAL  PROPERTY  BESCRIPTION 
REAB  INPUT  BATA 


WRITE (6 y 10) 

REAB ( 5  y / ) X I y  XMAX  y  BX  y  BT 
WRITE(6y20) 

REAB ( 5  y / ) SG I y  SGO  y  M 

PRINT  HEABING  ANB  INPUT  BATA 

WRITE ( 6  y  30 ) XI y  XMAX  y  SGO  y  SGI y  ETM  y  ERM  y  ATM  y  ARM  y  VTR  y  BT  y  M 
WRITE (6 y 40) 


